{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "private_outputs": true,
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/mikexcohen/Statistics_book/blob/main/stats_ch11_ttest.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Modern statistics: Intuition, Math, Python, R\n",
    "## Mike X Cohen (sincxpress.com)\n",
    "#### https://www.amazon.com/dp/B0CQRGWGLY\n",
    "#### Code for chapter 11 (t-tests)\n",
    "\n",
    "---\n",
    "\n",
    "# About this code file:\n",
    "\n",
    "### This notebook will reproduce most of the figures in this chapter (some figures were made in Inkscape), and illustrate the statistical concepts explained in the text. The point of providing the code is not just for you to recreate the figures, but for you to modify, adapt, explore, and experiment with the code.\n",
    "\n",
    "### Solutions to all exercises are at the bottom of the notebook.\n",
    "\n",
    "#### This code was written in google-colab. The notebook may require some modifications if you use a different IDE."
   ],
   "metadata": {
    "id": "mzfXc9E3Xq7k"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dFT1qwVEyxTW"
   },
   "outputs": [],
   "source": [
    "# import libraries and define global settings\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# pandas/seaborn for ex12+\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "# define global figure properties used for publication\n",
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats('svg') # display figures in vector format\n",
    "plt.rcParams.update({'font.size':14,             # font size\n",
    "                     'savefig.dpi':300,          # output resolution\n",
    "                     'axes.titlelocation':'left',# title location\n",
    "                     'axes.spines.right':False,  # remove axis bounding box\n",
    "                     'axes.spines.top':False,    # remove axis bounding box\n",
    "                     })"
   ]
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "BQlwO5QkJv7v"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.1: Goals of the t-test"
   ],
   "metadata": {
    "id": "YgKc7uVRJv5D"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "_,axs = plt.subplots(1,3,figsize=(10,3))\n",
    "\n",
    "## panel A: one-sample t-test\n",
    "data = np.random.normal(loc=.5,size=30)\n",
    "axs[0].plot(data,'ko',markersize=9,markerfacecolor=(.8,.8,.8))\n",
    "axs[0].plot([0,len(data)],[0,0],'k--',zorder=-10)\n",
    "axs[0].set(xlabel='Data index',ylabel='Data value',yticks=[],\n",
    "           title=r'$\\bf{A}$)  One sample')\n",
    "\n",
    "\n",
    "\n",
    "## panel B: paired-samples t-test\n",
    "N = 20\n",
    "data1 = np.random.normal(size=N)\n",
    "data2 = data1 + .5 + np.random.randn(N)*.4\n",
    "for x,y in zip(data1,data2):\n",
    "  # pick a random color\n",
    "  c = np.random.uniform(low=0,high=.8)\n",
    "\n",
    "  # plot it\n",
    "  axs[1].plot([0,1],[x,y],'o-',markersize=9,markerfacecolor=(c,c,c),color=(c,c,c))\n",
    "\n",
    "# plot adjustments\n",
    "axs[1].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=['pre','post'],\n",
    "           yticks=[],ylabel='Data value',title=r'$\\bf{B}$)  Paired samples')\n",
    "\n",
    "\n",
    "\n",
    "## panel C: two-samples t-test\n",
    "for i in range(2):\n",
    "\n",
    "  # create the data\n",
    "  data = np.random.normal(loc=i,scale=(i+1)/2,size=1000)\n",
    "\n",
    "  # take their histogram\n",
    "  yy,xx = np.histogram(data,bins='fd')\n",
    "  xx = (xx[1:]+xx[:-1])/2\n",
    "\n",
    "  # plot\n",
    "  c = i/2\n",
    "  axs[2].plot(xx,yy,linewidth=3,color=(c,c,c))\n",
    "\n",
    "# plot adjustments\n",
    "axs[2].set(xlabel='Exam score',ylabel='Count',xticks=[],yticks=[],title=r'$\\bf{C}$)  Two ind. samples')\n",
    "axs[2].legend(['Group 1','Group 2'],loc='upper right',fontsize=10)\n",
    "\n",
    "\n",
    "# final figure adjustments\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ttestGoals.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "M4_sU0UxJxuI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "YipkIzGLJvxV"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.2: A t-pdf"
   ],
   "metadata": {
    "id": "xvApsLWBosG0"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "t = np.linspace(-4,4,573)\n",
    "\n",
    "# a pdf with df=20\n",
    "tpdf = stats.t.pdf(t,20)\n",
    "\n",
    "plt.figure(figsize=(3,4))\n",
    "plt.plot(t,tpdf,'k',linewidth=2)\n",
    "plt.xlabel('T value')\n",
    "plt.ylabel('Probability')\n",
    "plt.yticks([])\n",
    "plt.ylim([-.01,np.max(tpdf)*1.05])\n",
    "plt.xlim(t[[0,-1]])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_tpdf.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "ezZ6iPcKovK7"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "ZeVa0JzikJ2d"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Computing p-values for one-tailed and two-tailed tests"
   ],
   "metadata": {
    "id": "3IdqtCoLktpF"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "tval = 2.1\n",
    "df = 13\n",
    "\n",
    "pvalL = stats.t.cdf(-tval,df)\n",
    "pvalR = stats.t.sf(tval,df)#1-stats.t.cdf(tval,df)\n",
    "pval2 = pvalR+pvalL\n",
    "\n",
    "print(f'One-tailed p-value on the left:  {pvalL}')\n",
    "print(f'One-tailed p-value on the right: {pvalR}')\n",
    "print(' ')\n",
    "print(f'Two-tailed p-value as the sum:   {pvalL+pvalR}')\n",
    "print(f'Two-tailed p-value by doubling:  {2*pvalL}')"
   ],
   "metadata": {
    "id": "-PHLJ-aFkJ0T"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# 1-cdf vs survival function:\n",
    "pvalC = 1-stats.t.cdf(tval,df)\n",
    "pvalS = stats.t.sf(tval,df) # sf = survival function\n",
    "\n",
    "print(f'P-value from 1-cdf: {pvalC}')\n",
    "print(f'P-value from s.f.:  {pvalS}')\n",
    "print(f'Difference:         {pvalC-pvalS}')\n",
    "\n",
    "# Conclusion: The difference for this particular t-value is at machine precision.\n",
    "#  Still, there's no harm in being slightly more accurate, so you can use sf instead of 1-cdf."
   ],
   "metadata": {
    "id": "GeTdpMKwtEwg"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "oIOAb8hRkJxQ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.3: T-values from p-values"
   ],
   "metadata": {
    "id": "X7lwZSnHris1"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# parameters\n",
    "t = np.linspace(-4,4,75)\n",
    "df = 13\n",
    "\n",
    "# cdf based on t-values\n",
    "cdf = stats.t.cdf(t,df)\n",
    "\n",
    "# t-values based on cdf\n",
    "pvals = np.linspace(.001,.999,73)\n",
    "tVals = -stats.t.isf(pvals,df) # sf is 1-cdf, and isf is the inverse of the sf. so 1-isf is inv(cdf)\n",
    "\n",
    "# note: the above line is the same as below, but isf has slightly higher accuracy\n",
    "#tVals = stats.t.ppf(pvals,df) # ppf = Percent point function, this is the inverse of the cdf\n",
    "\n",
    "\n",
    "# visualize\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "axs[0].plot(t,cdf,'k',linewidth=2)\n",
    "axs[0].set(xlabel='t value',ylabel='cdf',title=r'$\\bf{A}$)  CDF from t-values')\n",
    "\n",
    "axs[1].plot(pvals,tVals,'k',linewidth=2)\n",
    "axs[1].set(ylabel='t value',xlabel='cdf',title=r'$\\bf{B}$)  T-values from CDF')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_tFromP.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "lnzHK1o9rl9E"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# example usage to get the t-value associated with p=.05 and df=13\n",
    "pval = .05\n",
    "tFromP_L  = -stats.t.isf(  pval,df) # negative of the opposite-cdf for the left tail\n",
    "tFromP_R1 = -stats.t.isf(1-pval,df) # negative of the opposite-cdf for the right tail\n",
    "tFromP_R2 =  stats.t.isf(  pval,df) # opposite-cdf for the left tail\n",
    "\n",
    "print(f'Variable tFromP_L:  {tFromP_L:.3f}')\n",
    "print(f'Variable tFromP_R1: {tFromP_R1:.3f}')\n",
    "print(f'Variable tFromP_R2: {tFromP_R2:.3f}')"
   ],
   "metadata": {
    "id": "ZZI4QGCxriqV"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "uuifTy6uovNY"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.4: Example t-value"
   ],
   "metadata": {
    "id": "ph6ECiTwosD0"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# empirical t-value and df\n",
    "tval = 1.6\n",
    "df   = 20\n",
    "alpha = .05\n",
    "\n",
    "# redefine the t-values and corresponding pdf\n",
    "t = np.linspace(-4,4,573)\n",
    "tpdf = stats.t.pdf(t,20)\n",
    "\n",
    "\n",
    "# its associated p-value (but this is one-tailed for visualization; see text and next cell!)\n",
    "pval = 1-stats.t.cdf(tval,df)\n",
    "\n",
    "# critical t-value for alpha\n",
    "tCrit = stats.t.isf(alpha/2,df) # /2 for two-tailed!\n",
    "pHalf = np.max(tpdf)/2 # 1/2 max. (vertical) p(t), used for plotting\n",
    "\n",
    "\n",
    "\n",
    "plt.figure(figsize=(7,4))\n",
    "\n",
    "# plot the t distribution\n",
    "plt.plot(t,tpdf,'k',linewidth=1,label=r'$t_{20}$-pdf under H$_0$')\n",
    "\n",
    "# plot the dashed line for the critical t-value\n",
    "plt.axvline(tCrit,linestyle='--',color='gray')\n",
    "plt.text(tCrit-.02,pHalf*2,r'$\\alpha/2$ = %g'%(alpha/2),rotation=90,va='top',ha='right')\n",
    "\n",
    "# arrow and formula for the empirical t-value\n",
    "plt.gca().annotate(r'$t_{df}=\\frac{\\overline{x}-h_0}{s/\\sqrt{n}}$=%g'%tval,xytext=(tval+1,pHalf),\n",
    "                xy=(tval,0), xycoords='data',size=18,\n",
    "                arrowprops=dict(arrowstyle='->', color='k',linewidth=2,\n",
    "                                connectionstyle='angle,angleA=0,angleB=-90,rad=0'))\n",
    "\n",
    "\n",
    "# shaded area to the right of the empirical t-value\n",
    "tidx = np.argmin(np.abs(t-tval))\n",
    "plt.gca().fill_between(t[tidx:],tpdf[tidx:],color='k',alpha=.4)\n",
    "\n",
    "# and its annotation\n",
    "tidx = np.argmin(np.abs(t-(tval+t[-1])/2))\n",
    "plt.gca().annotate(f'{100*pval:.2f}%',xy=(t[tidx],tpdf[tidx]),\n",
    "            xytext=(t[tidx]+1,pHalf/2),ha='center',arrowprops={'color':'k'})\n",
    "\n",
    "\n",
    "# some final adjustments\n",
    "plt.xlabel('T value')\n",
    "plt.ylabel(r'$p(t|H_0)$')\n",
    "plt.yticks([])\n",
    "plt.ylim([-.01,np.max(tpdf)*1.05])\n",
    "plt.xlim([-1,t[-1]])\n",
    "plt.ylim([0,pHalf*2.1])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_tEmpWEq.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "qPfrhs6-p-rp"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "eJYHkffTp-un"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.5: Completion of the previous figure to show both tails"
   ],
   "metadata": {
    "id": "ApKXnHPA3Ed4"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(7,4))\n",
    "\n",
    "# plot the t distribution\n",
    "plt.plot(t,tpdf,'k',linewidth=1,label=r'$t_{20}$-pdf under H$_0$')\n",
    "\n",
    "# plot the dashed line for the critical t-value on the right side\n",
    "plt.axvline(tCrit,linestyle='--',color='gray')\n",
    "plt.text(tCrit-.02,pHalf*2,r'$\\alpha/2$ = %g'%(alpha/2),rotation=90,va='top',ha='right')\n",
    "\n",
    "# and again for the left side\n",
    "plt.axvline(-tCrit,linestyle='--',color='gray')\n",
    "plt.text(-tCrit+.028,pHalf*2,r'$\\alpha/2$ = %g'%(alpha/2),rotation=90,va='top',ha='left')\n",
    "\n",
    "# arrow and formula for the empirical t-value\n",
    "plt.gca().annotate(r'$t=$%g'%tval,xytext=(tval,pHalf),\n",
    "                xy=(tval,0), xycoords='data',size=18,ha='center',bbox=dict(fc='w',edgecolor='none'),\n",
    "                arrowprops=dict(arrowstyle='->', color='k',linewidth=2))\n",
    "\n",
    "# repeat on the left\n",
    "plt.gca().annotate(r'$t=$-%g'%tval,xytext=(-tval,pHalf),\n",
    "                xy=(-tval,0), xycoords='data',size=18,ha='center',bbox=dict(fc='w',edgecolor='none'),\n",
    "                arrowprops=dict(arrowstyle='->', color='k',linewidth=2))\n",
    "\n",
    "\n",
    "\n",
    "# shaded area to the right of the empirical t-value\n",
    "tidx = np.argmin(np.abs(t-tval))\n",
    "plt.gca().fill_between(t[tidx:],tpdf[tidx:],color='k',alpha=.4)\n",
    "tidx = np.argmin(np.abs(t--tval))\n",
    "plt.gca().fill_between(t[:tidx],tpdf[:tidx],color='k',alpha=.4)\n",
    "\n",
    "# and its annotation for the right side\n",
    "tidx = np.argmin(np.abs(t-(tval+t[-1])/2))\n",
    "plt.gca().annotate(f'{100*pval:.2f}%',xy=(t[tidx],tpdf[tidx]),xytext=(t[tidx]+1,pHalf/2),ha='center',arrowprops={'color':'k'})\n",
    "\n",
    "# now for the left side\n",
    "tidx = np.argmin(np.abs(t-(-tval+t[0])/2))\n",
    "plt.gca().annotate(f'{100*pval:.2f}%',xy=(t[tidx],tpdf[tidx]),xytext=(t[tidx]-.5,pHalf/2),ha='center',arrowprops={'color':'k'})\n",
    "\n",
    "\n",
    "# some final adjustments\n",
    "plt.xlabel('T value')\n",
    "plt.ylabel(r'$p(t|H_0)$')\n",
    "plt.yticks([])\n",
    "plt.ylim([-.01,np.max(tpdf)*1.05])\n",
    "plt.xlim(t[[0,-1]])\n",
    "plt.ylim([0,pHalf*2.1])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_tEmpWEq2.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "gdGDnXQNp-zo"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "uS2WorJ5osBE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.6: Testing for normality"
   ],
   "metadata": {
    "id": "MWJQ0jDy10Lt"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# the data\n",
    "data1 = np.random.randn(100)\n",
    "data2 = np.exp( np.random.randn(100) )\n",
    "\n",
    "# omnibus test\n",
    "Otest1 = stats.normaltest(data1)\n",
    "Otest2 = stats.normaltest(data2)\n",
    "\n",
    "# Shapiro's test\n",
    "Stest1 = stats.shapiro(data1)\n",
    "Stest2 = stats.shapiro(data2)\n",
    "\n",
    "# report the results\n",
    "print(f'Omnibus test in X1 (H0=normal): p={Otest1.pvalue:.3f}')\n",
    "print(f'Omnibus test in X2 (H0=normal): p={Otest2.pvalue:.3f}')\n",
    "print('')\n",
    "print(f'Shapiro test in X1 (H0=normal): p={Stest1.pvalue:.3f}')\n",
    "print(f'Shapiro test in X2 (H0=normal): p={Stest2.pvalue:.3f}')\n",
    "\n",
    "# show the histograms\n",
    "yy1,xx1 = np.histogram(data1,bins='fd')\n",
    "xx1 = (xx1[1:]+xx1[:-1])/2\n",
    "yy2,xx2 = np.histogram(data2,bins='fd')\n",
    "xx2 = (xx2[1:]+xx2[:-1])/2\n",
    "\n",
    "\n",
    "# plotting\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.plot(xx1,yy1,'k--',linewidth=3,label=r'$X_1$')\n",
    "plt.plot(xx2,yy2,linewidth=3,color=(.5,.5,.5),label=r'$X_2$')\n",
    "plt.gca().set(xlabel='Data value',ylabel='Count')\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_normTests.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "isvJapnV111V"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "AKW3W9AF10DL"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.7: Increasing the t-value"
   ],
   "metadata": {
    "id": "ASmY9pVxor-H"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "x = np.linspace(-4,4,501)\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,4))\n",
    "\n",
    "\n",
    "## panel A: probably not significant\n",
    "g1 = stats.norm.pdf(x,loc=-.3,scale=1)\n",
    "g2 = stats.norm.pdf(x,loc= .3,scale=1)\n",
    "\n",
    "axs[0].plot(x,g1,color='k',linewidth=2)\n",
    "axs[0].plot(x,g2,'--',color=(.6,.6,.6),linewidth=2)\n",
    "axs[0].set(xticks=[],xlim=x[[0,-1]],yticks=[],ylabel='Probability',\n",
    "           title=r'$\\bf{A}$)  Non-significant')\n",
    "\n",
    "\n",
    "## panel B: significant by larger mean difference\n",
    "g1 = stats.norm.pdf(x,loc=-1,scale=1)\n",
    "g2 = stats.norm.pdf(x,loc= 1,scale=1)\n",
    "\n",
    "axs[1].plot(x,g1,color='k',linewidth=2)\n",
    "axs[1].plot(x,g2,'--',color=(.6,.6,.6),linewidth=2)\n",
    "axs[1].set(xticks=[],xlim=x[[0,-1]],yticks=[],ylabel='Probability',\n",
    "           title=r'$\\bf{B}$)  Large mean distance')\n",
    "\n",
    "\n",
    "## panel C: significant by reduced variance\n",
    "g1 = stats.norm.pdf(x,loc=-.3,scale=.2)\n",
    "g2 = stats.norm.pdf(x,loc= .3,scale=.2)\n",
    "\n",
    "axs[2].plot(x,g1,color='k',linewidth=2)\n",
    "axs[2].plot(x,g2,'--',color=(.6,.6,.6),linewidth=2)\n",
    "axs[2].set(xticks=[],xlim=x[[0,-1]],yticks=[],xlabel='Data value',ylabel='Probability',\n",
    "           title=r'$\\bf{C}$)  Low variance')\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_sigMecs.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "9VBwa_Ef3T4h"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "LqydXjgvQ2E3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# One-sample t-test"
   ],
   "metadata": {
    "id": "q7s93I3gQ2CI"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# given data\n",
    "X = np.array([80, 85, 90, 70, 75, 72, 88, 77, 82, 65, 79, 81, 74, 86, 68])\n",
    "h0 = 75\n",
    "\n",
    "# descriptives\n",
    "meanX = np.mean(X)\n",
    "stdX  = np.std(X,ddof=1)\n",
    "ssize = len(X)\n",
    "\n",
    "# t-value\n",
    "tval = (meanX-h0) / (stdX/np.sqrt(ssize))\n",
    "\n",
    "# p-value\n",
    "pval = 1-stats.t.cdf(tval,ssize-1)\n",
    "pval *= 2 # two-tailed!\n",
    "\n",
    "# print everything out!\n",
    "print(f'Sample mean: {meanX:.2f}')\n",
    "print(f'Sample std:  {stdX:.2f}')\n",
    "print(f'Sample size: {ssize}')\n",
    "print('')\n",
    "print(f'T-value: {tval:.3f}')\n",
    "print(f'p-value: {pval:.3f}')"
   ],
   "metadata": {
    "id": "yr0Ql-uhQ1_b"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Repeat using the stats libary:\n",
    "ttest = stats.ttest_1samp(X,h0)\n",
    "\n",
    "# the output variable is its own type\n",
    "print( type(ttest) )\n",
    "\n",
    "# which contains three elements:\n",
    "print('')\n",
    "print(ttest)\n",
    "\n",
    "# let's print the results\n",
    "print('')\n",
    "print('Results from stats.ttest_1samp:')\n",
    "print(f't({ttest.df})={ttest.statistic:.3f}, p<{ttest.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "hLnN0NmzQ18i"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# btw, data are consistent with a normal distribution\n",
    "print(f'Shapiro p-value = {stats.shapiro(X).pvalue:.2f}')"
   ],
   "metadata": {
    "id": "ieNAgdXnQ13F"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "Lco2nRlG60_o"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.8: Paired-samples t-test"
   ],
   "metadata": {
    "id": "lvOrmWkqQ10U"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# the data\n",
    "Xn = np.array([ 60, 52, 90, 20, 33, 95, 18, 47, 78, 65 ])\n",
    "Xq = np.array([ 65, 60, 84, 23, 37, 95, 17, 53, 88, 66 ])\n",
    "sampsize = len(Xn)\n",
    "\n",
    "# their difference\n",
    "Delta = Xq-Xn\n",
    "\n",
    "# visualize\n",
    "_,axs = plt.subplots(1,2,figsize=(7,3))\n",
    "\n",
    "## draw the individual lines\n",
    "for i,j in zip(Xn,Xq):\n",
    "  axs[0].plot([0,1],[i,j],'o-',color=(.8,.8,.8),\n",
    "              markersize=12,markerfacecolor=(.8,.8,.8),markeredgecolor='k')\n",
    "\n",
    "axs[0].set(xlim=[-1,2],xticks=[0,1],ylabel='Scores',xticklabels=[r'X$_N$',r'X$_Q$'],\n",
    "           ylim=[0,100],title=r'$\\bf{A}$)  Raw data')\n",
    "\n",
    "\n",
    "## draw the difference scores\n",
    "axs[1].plot(np.zeros(sampsize),Delta,'ko',markersize=12,markerfacecolor=(.8,.8,.8))\n",
    "axs[1].plot([-.1,.1],[0,0],'k--',zorder=-1)\n",
    "axs[1].set(xticks=[0],ylabel='Difference scores',xticklabels=[r'$\\Delta$'],\n",
    "           ylim=[-50,50],title=r'$\\bf{B}$)  Differences')\n",
    "\n",
    "# export\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_pairedTtest.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "3bXMQKChhqRD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# test it!\n",
    "ttest = stats.ttest_1samp(Delta,0)\n",
    "\n",
    "# print the results\n",
    "print(f't({ttest.df})={ttest.statistic:.3f}, p<{ttest.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "mhdDiKZohqT-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# btw, data are consistent with a normal distribution\n",
    "print(f'Xn Shapiro p-value = {stats.shapiro(Xn).pvalue:.2f}')\n",
    "print(f'Xq Shapiro p-value = {stats.shapiro(Xq).pvalue:.2f}')\n",
    "print(f'Xy Shapiro p-value = {stats.shapiro(Delta).pvalue:.2f}')"
   ],
   "metadata": {
    "id": "NaGUcyRvhqWg"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "-Uq4-Jdl70wT"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.9: Example of 2-sample ttest"
   ],
   "metadata": {
    "id": "PiMFTrrSQ1u3"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate data\n",
    "data1 = stats.exponnorm.rvs(3,size=50)\n",
    "data2 = stats.gumbel_r.rvs(size=42)\n",
    "\n",
    "# compute their histograms\n",
    "yy1,xx1 = np.histogram(data1,bins='fd')\n",
    "xx1 = (xx1[1:]+xx1[:-1])/2\n",
    "yy2,xx2 = np.histogram(data2,bins='fd')\n",
    "xx2 = (xx2[1:]+xx2[:-1])/2\n",
    "\n",
    "\n",
    "\n",
    "# show the data!\n",
    "_,axs = plt.subplots(1,2,figsize=(7,3.5))\n",
    "\n",
    "## raw data\n",
    "axs[0].plot(np.random.randn(len(data1))/40,data1,\n",
    "            'ko',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[0].plot(np.random.randn(len(data2))/40+1,data2,\n",
    "            'ks',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[0].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=[r'$X_1$',r'$X_2$'],\n",
    "           title=r'$\\bf{A}$)  Data')\n",
    "\n",
    "\n",
    "## histograms\n",
    "axs[1].plot(xx1,yy1,'ko--',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6,linewidth=2,label=r'$X_1$')\n",
    "axs[1].plot(xx2,yy2,'ks-',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6,linewidth=2,label=r'$X_2$')\n",
    "axs[1].set(xlabel='Data value',ylabel='Count',title=r'$\\bf{B}$)  Distributions')\n",
    "axs[1].legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_indTtest.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "EhLUhOX7pToI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# doubling rubric\n",
    "s1 = np.std(data1,ddof=1)\n",
    "s2 = np.std(data2,ddof=1)\n",
    "\n",
    "# report\n",
    "print(f'Standard deviations are {s1:.2f} and {s2:.2f}')\n",
    "print(f'Ratio of max:min stdevs is {np.max([s1,s2])/np.min([s1,s2]):.2f}')\n",
    "\n",
    "# Levene's test\n",
    "lres = stats.levene(data1,data2)\n",
    "print('')\n",
    "print(f\"Levene's test for homogeneity of variance: W={lres.statistic:.2f}, p={lres.pvalue:.3f}\")"
   ],
   "metadata": {
    "id": "gnjg72jZQ1sO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## tests for normal distribution\n",
    "\n",
    "# omnibus test\n",
    "Otest1 = stats.normaltest(data1)\n",
    "Otest2 = stats.normaltest(data2)\n",
    "\n",
    "print(f'Omnibus test in X1 (H0=normal): p={Otest1.pvalue:.3f}')\n",
    "print(f'Omnibus test in X2 (H0=normal): p={Otest2.pvalue:.3f}')\n",
    "print('')\n",
    "\n",
    "\n",
    "# Shapiro's test\n",
    "Stest1 = stats.shapiro(data1)\n",
    "Stest2 = stats.shapiro(data2)\n",
    "\n",
    "print(f'Shapiro test in X1 (H0=normal): p={Stest1.pvalue:.3f}')\n",
    "print(f'Shapiro test in X2 (H0=normal): p={Stest2.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "1TRzyA8Sw3Kv"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# now for the t-test\n",
    "tres = stats.ttest_ind(data1,data2,equal_var=False)\n",
    "print(f't={tres.statistic:.2f}, p={tres.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "EicajwSSuWIP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# FYI, here's the result assuming equal variance (see also Exercise 9)\n",
    "tres = stats.ttest_ind(data1,data2,equal_var=True)\n",
    "print(f't={tres.statistic:.2f}, p={tres.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "Op_k-tM8or7H"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "GVl4kD4XOtO7"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11:10: Wilcoxon signed-rank"
   ],
   "metadata": {
    "id": "ALNZfi53OtMe"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# the data\n",
    "data = np.random.randn(100)**2\n",
    "h0 = 1\n",
    "\n",
    "# show the data!\n",
    "_,axs = plt.subplots(1,2,figsize=(7,3.5))\n",
    "\n",
    "## raw data\n",
    "axs[0].plot(data,'ko',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[0].axhline(h0,linestyle='--',color=(.4,.4,.4),linewidth=2)\n",
    "axs[0].axhline(np.median(data),color='k',linewidth=3)\n",
    "axs[0].set(xlabel='Data index',ylabel='Data value',title=r'$\\bf{A}$)  Data')\n",
    "\n",
    "## histogram\n",
    "axs[1].hist(data,bins='fd',facecolor=(.9,.9,.9),edgecolor='k',label='Data')\n",
    "axs[1].axvline(h0,linestyle='--',color=(.4,.4,.4),linewidth=2,label=r'H$_0$')\n",
    "axs[1].axvline(np.median(data),color='k',linewidth=3,label='Median')\n",
    "axs[1].set(xlabel='Data value',ylabel='Count',title=r'$\\bf{B}$)  Distribution')\n",
    "axs[1].legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ranktest.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "rWPhszB-O1Sr"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# the test!\n",
    "wtest = stats.wilcoxon(data-h0,method='approx')\n",
    "\n",
    "# and print the results\n",
    "print(f'Wilcoxon test: z={wtest.zstatistic:.2f}, p={wtest.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "MgjuHwOlO1Vi"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "RMhBLKqQOtHK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 11.11: Margin figure about the sign of z"
   ],
   "metadata": {
    "id": "knegbKSX6iFP"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the data, shifted by H0=1\n",
    "_,axs = plt.subplots(2,1,figsize=(3,6))\n",
    "\n",
    "for i in range(2):\n",
    "\n",
    "  # create and shift data\n",
    "  d = np.random.randn(30)\n",
    "  d += i*2-1\n",
    "\n",
    "  # Wilcoxon z-score\n",
    "  z = stats.wilcoxon(d,method='approx',alternative='two-sided').zstatistic\n",
    "\n",
    "  # draw the figure\n",
    "  axs[i].plot(d,range(len(d)),'ko',markersize=8)\n",
    "  axs[i].axvline(0,zorder=-1,color='gray')\n",
    "  axs[i].set(xlabel='Data values',ylabel='Data index',yticks=[],xlim=[-3,3])\n",
    "  axs[i].set_title(f'Wilcoxon z={z:.2f}',loc='center')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_wilcoxonSign.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "dRWtoKqLvL42"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Mann-Whitney U test"
   ],
   "metadata": {
    "id": "4uqt2O7iOtEj"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# same data as we used for the independent-samples t-test\n",
    "data1 = stats.exponnorm.rvs(3,size=50)\n",
    "data2 = stats.gumbel_r.rvs(size=42)\n",
    "\n",
    "# MW-U test\n",
    "mwu = stats.mannwhitneyu(data1,data2)\n",
    "print(f'U = {mwu.statistic:.0f}, p = {mwu.pvalue:.3f}')\n",
    "\n",
    "\n",
    "# parametric t-test (gives the same statistical conclusion as the MWU)\n",
    "tres = stats.ttest_ind(data1,data2,equal_var=False)\n",
    "print(f't = {tres.statistic:.2f}, p = {tres.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "MzeGoUiPOtCB"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "iE1ZBZmg0ppX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "OHZT_OCXor4I"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# parameters\n",
    "N  = 50\n",
    "h0 = -np.pi/2\n",
    "\n",
    "# create the dataset\n",
    "X = stats.laplace_asymmetric.rvs(2,size=N)\n",
    "dataMean = np.mean(X)\n",
    "\n",
    "# visualize the data\n",
    "_,axs = plt.subplots(1,2,figsize=(9,3))\n",
    "\n",
    "axs[0].plot(X,'kp',markersize=8,markerfacecolor=(.9,.9,.9),label='Data')\n",
    "axs[0].plot([0,N],[h0,h0],'k--',zorder=-10,linewidth=3,label=r'H$_0$ value')\n",
    "axs[0].plot([0,N],[dataMean,dataMean],'k:',linewidth=3,label='Emp. mean')\n",
    "axs[0].set(xlabel='Data index',ylabel='Data value')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Raw data')\n",
    "\n",
    "axs[1].hist(X,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[1].axvline(h0,linestyle='--',color='k',linewidth=3,label=r'H$_0$ value')\n",
    "axs[1].axvline(dataMean,linestyle=':',color='k',linewidth=3,label=r'Emp. mean')\n",
    "axs[1].set(xlabel='Data value',ylabel='Count')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Histogram and means')\n",
    "axs[1].legend(bbox_to_anchor=[1,.9])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex1.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Yzex2bxPLYKl"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# now for the t-tests\n",
    "\n",
    "## manual calculation\n",
    "t_num = dataMean - h0\n",
    "t_den = np.std(X,ddof=1) / np.sqrt(N)\n",
    "\n",
    "tval  = t_num / t_den\n",
    "pval  = 1-stats.t.cdf( np.abs(tval) ,N-1)\n",
    "pval *= 2 # double it for 2-tailed test\n",
    "\n",
    "\n",
    "## using scipy.stats\n",
    "r  = stats.ttest_1samp(X,h0)\n",
    "t  = r.statistic\n",
    "df = r.df\n",
    "p  = r.pvalue\n",
    "\n",
    "\n",
    "# print both results\n",
    "print(f'Manual ttest: t({N-1})={tval:.3f}, p={pval:.3f}')\n",
    "print(f'Scipy  ttest: t({df})={t:.3f}, p={p:.3f}')"
   ],
   "metadata": {
    "id": "uPl2GaD7or1N"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "Wiv9J-hBVigA"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "39HpgAP7VibS"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# how often do we get subthreshold results?\n",
    "\n",
    "nExps = 500\n",
    "issig = np.zeros(nExps,dtype=bool) # variable type 'bool' for convenience in plotting\n",
    "means = np.zeros(nExps)\n",
    "stds  = np.zeros(nExps)\n",
    "\n",
    "# run the experiment\n",
    "#   (Note: For a small extra challenge, you could re-implement this without\n",
    "#          a for-loop using matrix input, after completing the next exercise.)\n",
    "for i in range(nExps):\n",
    "\n",
    "  # generate data and store the mean/std\n",
    "  X = stats.laplace_asymmetric.rvs(2,size=N)\n",
    "  means[i] = np.mean(X)\n",
    "  stds[i]  = np.std(X,ddof=1)\n",
    "\n",
    "  # run the ttest and store if \"significant\"\n",
    "  r = stats.ttest_1samp(X,h0)\n",
    "  issig[i] = r.pvalue<.05\n",
    "\n",
    "# print the results\n",
    "print(f'p<.05 in {np.sum(issig)}/{nExps} times.')\n",
    "\n",
    "\n",
    "\n",
    "# show the results\n",
    "_,axs = plt.subplots(1,2,figsize=(7,3))\n",
    "\n",
    "# means\n",
    "axs[0].plot(np.random.randn(sum(issig))/40,means[issig],\n",
    "            'ko',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[0].plot(np.random.randn(sum(~issig))/40+1,means[~issig],\n",
    "            'ks',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[0].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=['p<.05','p>.05'],\n",
    "           title=r'$\\bf{A}$)  Sample means')\n",
    "\n",
    "# stds\n",
    "axs[1].plot(np.random.randn(sum(issig))/40,stds[issig],\n",
    "            'ko',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[1].plot(np.random.randn(sum(~issig))/40+1,stds[~issig],\n",
    "            'ks',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.6)\n",
    "axs[1].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=['p<.05','p>.05'],\n",
    "           title=r'$\\bf{B}$)  Sample stds.')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex2.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "CiwrS7ymoryP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "_veJ-lcaVfvX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 3"
   ],
   "metadata": {
    "id": "lUHoAM4BVfsX"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "NperSample = 40\n",
    "MDatasets = 25\n",
    "\n",
    "# data\n",
    "X = np.random.normal(loc=1,scale=1,size=(NperSample,MDatasets))\n",
    "\n",
    "# confirm data size\n",
    "print('Data size should be sample-size X datasets:')\n",
    "print(X.shape)"
   ],
   "metadata": {
    "id": "JpIMFAixX5TS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# ttest with matrix input\n",
    "ttest_matrix = stats.ttest_1samp(X,0)\n",
    "\n",
    "# ttest in for-loop over each column (each dataset)\n",
    "ttest_4loop = np.zeros(MDatasets)\n",
    "for i in range(MDatasets):\n",
    "  ttest_4loop[i] = stats.ttest_1samp(X[:,i],0).statistic\n",
    "\n",
    "\n",
    "# print the results\n",
    "print('Matrix  |  Vector')\n",
    "print('--------|--------')\n",
    "for i in range(MDatasets):\n",
    "  print(f'{ttest_matrix.statistic[i]:.4f}  |  {ttest_4loop[i]:.4f}')"
   ],
   "metadata": {
    "id": "rbc7FGbpX5Vc"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "xYNzHtImkwIx"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "VhLKcFCbkwF_"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# data parameters\n",
    "N = 40\n",
    "k = 300\n",
    "\n",
    "# list of standard deviations\n",
    "stds = np.linspace(.1,3,k)\n",
    "\n",
    "# initialize the t/p vectors\n",
    "t = np.zeros(k)\n",
    "p = np.zeros(k)\n",
    "s = np.zeros(k) # this line is for exercise 5\n",
    "\n",
    "for i in range(len(stds)):\n",
    "  X = np.random.normal(0,stds[i],size=N)\n",
    "  X = X-np.mean(X) + .5 # force mean=.5\n",
    "  ttest = stats.ttest_1samp(X,0)\n",
    "  t[i]  = ttest.statistic\n",
    "  p[i]  = ttest.pvalue\n",
    "\n",
    "  # get the sample std (used in exercise 5)\n",
    "  s[i] = np.std(X,ddof=1)\n",
    "\n",
    "\n",
    "# and now the plotting\n",
    "_,axs = plt.subplots(1,3,figsize=(10,3))\n",
    "\n",
    "# t's\n",
    "axs[0].plot(stds,t,'ks',markerfacecolor='w')\n",
    "axs[0].set(xlabel='Standard deviation',ylabel='t-value',title=r'$\\bf{A}$)  T-values')\n",
    "\n",
    "# p's\n",
    "axs[1].plot(stds,p,'ks',markerfacecolor='w')\n",
    "axs[1].set(xlabel='Standard deviation',ylabel='p-value',title=r'$\\bf{B}$)  P-values')\n",
    "\n",
    "# t and p\n",
    "axs[2].plot(t,p,'ks',markerfacecolor='w')\n",
    "axs[2].set(xlabel='T-value',ylabel='p-value',title=r'$\\bf{C}$)  p by t')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex4.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "BZKbZt3LkxdG"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "qPXn8Q9fBDm0"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "JWSthVKeBDkY"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# No, it doesn't really matter, because even with N=40, the sample standard deviation\n",
    "# is a fairly accurate estimate of the population standard deviation, certainly for\n",
    "# this range of standard deviation values.\n",
    "\n",
    "# You can recreate the figure by replacing variable 'stds' with 's' in the code above,\n",
    "# and by demonstrating the strong correlation between sample and theoretical standard deviation.\n",
    "\n",
    "# correlation coefficient (values close to 1 indicate a very strong relationship)\n",
    "r = np.corrcoef(stds,s)\n",
    "\n",
    "# plot\n",
    "plt.plot(stds,s,'ko')\n",
    "plt.xlabel('Theoretical population standard deviations')\n",
    "plt.ylabel('Empirical sample standard deviations')\n",
    "plt.title(f'Correlation: r={r[0,1]:.3f}',loc='center')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "EWdb3p-Fm4dA"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "_6ewsNgRB2aI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "DXxbpUU6qmQY"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "nExperiments = 250\n",
    "meanoffsets = np.linspace(0,.3,51)\n",
    "samplesizes = np.arange(10,811,step=50)\n",
    "\n",
    "\n",
    "# initialize\n",
    "propSig = np.zeros((len(samplesizes),len(meanoffsets)))\n",
    "\n",
    "\n",
    "# loop over sample sizes\n",
    "for sidx,ssize in enumerate(samplesizes):\n",
    "\n",
    "  # loop over effect sizes\n",
    "  for eidx,effect in enumerate(meanoffsets):\n",
    "\n",
    "    # generate the data\n",
    "    X = np.random.normal(loc=effect,scale=1.5,size=(ssize,nExperiments))\n",
    "\n",
    "    # run the t-test and store the results\n",
    "    T = stats.ttest_1samp(X,0)\n",
    "    propSig[sidx,eidx] = 100*np.mean( T.pvalue<.05 )\n",
    "\n",
    "\n",
    "# visualize in a matrix\n",
    "plt.imshow(propSig,extent=[meanoffsets[0],meanoffsets[-1],samplesizes[0],samplesizes[-1]],\n",
    "           vmin=0,vmax=25,origin='lower',aspect='auto',cmap='gray')\n",
    "plt.xlabel('Mean offset')\n",
    "plt.ylabel('Sample size')\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('Percent tests with p<.05')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex6.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "9YQpQl4vm4aK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "QnntaIMEVfpp"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "gWnMCZbsVfmz"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "Xn = np.array([ 60, 52, 90, 20, 33, 95, 18, 47, 78, 65 ])\n",
    "Xq = np.array([ 65, 60, 84, 23, 37, 95, 17, 53, 88, 66 ])\n",
    "sampsize = len(Xn)\n",
    "\n",
    "# simple subtraction (Y1 in the text)\n",
    "Ysub = Xq-Xn\n",
    "\n",
    "# zscore subtraction (Y2 in the text)\n",
    "Ysbz = stats.zscore(Xq) - stats.zscore(Xn)\n",
    "\n",
    "# percent change (Y3 in the text)\n",
    "Ypct = 100*(Xq-Xn) / Xn\n",
    "\n",
    "# normalized ratio (Y4 in the text)\n",
    "Ynrt = (Xq-Xn) / (Xq+Xn)\n",
    "\n",
    "\n",
    "# plot\n",
    "_,axs = plt.subplots(2,3,figsize=(10,6))\n",
    "axs[0,0].plot(Ysub,Ysbz,'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[0,0].set(xlabel='Subtraction',ylabel='Z-scored')\n",
    "\n",
    "axs[0,1].plot(Ysub,Ypct,'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[0,1].set(xlabel='Subtraction',ylabel='Percent change')\n",
    "\n",
    "axs[0,2].plot(Ysub,Ynrt,'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[0,2].set(xlabel='Subtraction',ylabel='Norm. ratio')\n",
    "\n",
    "axs[1,0].plot(Ysbz,Ypct,'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[1,0].set(xlabel='Z-scored',ylabel='Percent change')\n",
    "\n",
    "axs[1,1].plot(Ysbz,Ynrt,'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[1,1].set(xlabel='Z-scored',ylabel='Norm. ratio')\n",
    "\n",
    "axs[1,2].plot(Ypct,Ynrt,'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[1,2].set(xlabel='Percent change',ylabel='Norm. ratio')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex7.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "e-Xb8vJiOQGx"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# t-tests\n",
    "tSub = stats.ttest_1samp(Ysub,0)\n",
    "tPct = stats.ttest_1samp(Ypct,0)\n",
    "tsbz = stats.ttest_1samp(Ysbz,0)\n",
    "tnrt = stats.ttest_1samp(Ynrt,0)\n",
    "\n",
    "# print the results\n",
    "print(f'Subtraction (Y1): t({tSub.df})={tSub.statistic:.3f}, p<{tSub.pvalue:.3f}')\n",
    "print(f'Percent chg (Y2): t({tPct.df})={tPct.statistic:.3f}, p<{tPct.pvalue:.3f}')\n",
    "print(f'Z subtract  (Y3): t({tsbz.df})={tsbz.statistic:.3f}, p<{tsbz.pvalue:.3f}')\n",
    "print(f'Norm. ratio (Y4): t({tnrt.df})={tnrt.statistic:.3f}, p<{tnrt.pvalue:.3f}')"
   ],
   "metadata": {
    "id": "KJgczvVUOP-a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "tLdNsSjgOP71"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "SSXp2TJ6OP5g"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# parameters\n",
    "mu1 = 1.2 # population mean in dataset 1\n",
    "mu2 = 1   # population mean in dataset 2\n",
    "\n",
    "# sample sizes\n",
    "ns = np.arange(10,201,step=10)\n",
    "\n",
    "\n",
    "# setup the figure\n",
    "_,axs = plt.subplots(2,1,figsize=(8,6))\n",
    "\n",
    "# start the experiment!\n",
    "for ni,N in enumerate(ns):\n",
    "\n",
    "  # generate the data (100 datasets at a time)\n",
    "  data1 = np.random.normal(loc=mu1,scale=.5,size=(N,100))\n",
    "  data2 = np.random.normal(loc=mu2,scale=.5,size=(N,100))\n",
    "\n",
    "  # run the ttest\n",
    "  ttest = stats.ttest_ind(data1,data2)\n",
    "  t = ttest.statistic;\n",
    "  p = ttest.pvalue;\n",
    "\n",
    "  # plot the t-value, colored by significance\n",
    "  axs[0].plot(np.full(np.sum(p>.05),N),t[p>.05],'ks',markersize=8,markerfacecolor=(.5,.5,.5),alpha=.3)\n",
    "  axs[0].plot(np.full(np.sum(p<.05),N),t[p<.05],'ro',markersize=8,markerfacecolor=(.7,.3,.3))\n",
    "\n",
    "  # plot the p-values\n",
    "  axs[1].plot(np.full(np.sum(p>.05),N),p[p>.05],'ks',markersize=8,markerfacecolor=(.5,.5,.5),alpha=.3)\n",
    "  axs[1].plot(np.full(np.sum(p<.05),N),p[p<.05],'ro',markersize=8,markerfacecolor=(.7,.3,.3))\n",
    "\n",
    "\n",
    "\n",
    "## rest of the visualization\n",
    "axs[0].set(xlabel='Sample size (equal groups)',xticks=ns[::2],ylabel='T-test value')\n",
    "axs[0].set_title(r'$\\bf{A}$)  T-values')\n",
    "\n",
    "# adjust the p-values panel\n",
    "axs[1].set(xlabel='Sample size (equal groups)',xticks=ns[::2],ylabel='P value')\n",
    "axs[1].set_yscale('log')\n",
    "axs[1].axhline(.05,linestyle='--',color='r')\n",
    "axs[1].set_title(r'$\\bf{B}$)  P-values')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex8.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "hhQVCop9VfkB"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# compute critical t-values for the degrees of freedom\n",
    "tCrit = stats.t.isf(.05/2,2*ns-2)\n",
    "\n",
    "# and visualize\n",
    "plt.figure(figsize=(7,3))\n",
    "plt.plot(ns,tCrit,'ko',markersize=10,markerfacecolor=(.6,.6,.6))\n",
    "plt.ylim([1.8,2.2])\n",
    "plt.xlabel('Degrees of freedom')\n",
    "plt.ylabel('Critical t-values (2-tailed)')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "yqOXoLzk3RLo"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "v9rWgB6zVfeU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 9"
   ],
   "metadata": {
    "id": "Wgz8JaA9VfbG"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# range of standard deviations\n",
    "stdevs = np.linspace(.01,15,41)\n",
    "\n",
    "# initialize results matrix\n",
    "results = np.zeros((3,len(stdevs)))\n",
    "tCrit = np.zeros(len(stdevs))\n",
    "\n",
    "# start the experiment!\n",
    "for si,std in enumerate(stdevs):\n",
    "\n",
    "  # create two groups of data\n",
    "  X1 = np.random.normal(loc=1,scale=1,size=50)\n",
    "  X2 = np.random.normal(loc=1.1,scale=std,size=40)\n",
    "\n",
    "  # levene's test\n",
    "  results[0,si] = np.log( stats.levene(X1,X2).pvalue )\n",
    "\n",
    "  # difference of t-values\n",
    "  same_var = stats.ttest_ind(X1,X2,equal_var=True)  # equal variance\n",
    "  diff_var = stats.ttest_ind(X1,X2,equal_var=False) # unequal variance\n",
    "  results[1,si] = same_var.statistic # equal variance\n",
    "  results[2,si] = diff_var.statistic # unequal variance\n",
    "\n",
    "\n",
    "  # compute df for tCrit\n",
    "  s1,s2 = np.var(X1,ddof=1),np.var(X2,ddof=1)\n",
    "  n1,n2 = len(X1),len(X2)\n",
    "  df_num = (s1/n1 + s2/n2)**2\n",
    "  df_den = s1**2/(n1**2*(n1-1)) + s2**2/(n2**2*(n2-1))\n",
    "\n",
    "  tCrit[si] = stats.t.isf(.05/2,df_num/df_den)\n",
    "\n",
    "\n",
    "\n",
    "# plot\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "# levene's test results\n",
    "axs[0].plot(stdevs,results[0,:],'ks',markersize=10,markerfacecolor='gray')\n",
    "axs[0].axhline(np.log(.05),color=(.6,.6,.6),linestyle='--',zorder=-1)\n",
    "axs[0].text(np.max(stdevs),np.log(.1),'p=.05',ha='right',color=(.6,.6,.6))\n",
    "axs[0].set(xlabel=r'Standard deviation of $X_2$',ylabel='log(p)',title=r\"$\\bf{A}$)  Levene's P-values\")\n",
    "\n",
    "# t-tests\n",
    "axs[1].plot(stdevs,results[1,:],'ks',markersize=8,markerfacecolor=(.4,.4,.4),label='Equal var.')\n",
    "axs[1].plot(stdevs,results[2,:],'ko',markersize=8,markerfacecolor=(.8,.8,.8),label='Unequal var.')\n",
    "axs[1].plot(stdevs,tCrit,'--',color=(.6,.6,.6),zorder=-1,label='p=.05')\n",
    "axs[1].plot(stdevs,-tCrit,'--',color=(.6,.6,.6),zorder=-1)\n",
    "axs[1].set(xlabel=r'Standard deviation of $X_2$',ylabel='T-value',title=r\"$\\bf{B}$)  T-test t-values\")\n",
    "axs[1].legend(fontsize=10,bbox_to_anchor=[.8,1])\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex9.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "cvPyqj-zujIR"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Not in the instructions, but I think it's also interesting to\n",
    "#  plot the ratio of t-values as a function of standard deviation\n",
    "# values >1 indicate a larger t-value for equal compared to unequal variance formula\n",
    "plt.plot(stdevs,results[1,:]/results[2,:],'ko',markersize=8)\n",
    "plt.axhline(1,linestyle='--',color='gray',zorder=-1)\n",
    "plt.ylim([.5,1.5])\n",
    "plt.ylabel('Equal to Unequal variance t ratio')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "CgE0DIpGvMZh"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# If you're bored on a Saturday night, you can also use this code\n",
    "# to explore the impact of sample sizes and mean offsets."
   ],
   "metadata": {
    "id": "R5luU0U2HwTE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "Au6DO8lvujZD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 10"
   ],
   "metadata": {
    "id": "EsyuBZcCujb3"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate the data\n",
    "sigmas = np.linspace(.1,1.2,20)\n",
    "\n",
    "# null hypothesis value\n",
    "h0 = .5\n",
    "\n",
    "# initialize the results matrices\n",
    "tvals = np.zeros((2,len(sigmas)))\n",
    "cents = np.zeros((2,len(sigmas)))\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(11,3))\n",
    "\n",
    "\n",
    "\n",
    "# compute and store all moments in a matrix\n",
    "for i,s in enumerate(sigmas):\n",
    "\n",
    "  # generate mean-centered data\n",
    "  X = np.exp(np.random.randn(100)*s)\n",
    "  X -= np.mean(X)\n",
    "\n",
    "  # compute and store the descriptives\n",
    "  cents[0,i] = np.mean(X) - h0\n",
    "  cents[1,i] = np.median(X) - h0\n",
    "\n",
    "  # draw the histogram\n",
    "  if i%3==0:\n",
    "    mc = len(sigmas)+2\n",
    "    y,x = np.histogram(X,bins='fd')\n",
    "    axs[0].plot((x[:-1]+x[1:])/2,y,color=(i/mc,i/mc,i/mc),linewidth=2)\n",
    "    axs[0].axvline(np.median(X),color=(i/mc,i/mc,i/mc),linestyle='--',linewidth=.8)\n",
    "\n",
    "  # parametric t-test\n",
    "  tvals[0,i] = stats.ttest_1samp(X,h0).statistic\n",
    "\n",
    "  # Wilcoxon test\n",
    "  tvals[1,i] = stats.wilcoxon(X-h0,method='approx').zstatistic\n",
    "\n",
    "\n",
    "## plot\n",
    "axs[0].set(xlim=[-1.5,4],xlabel='Data value',ylabel='Count')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Distributions')\n",
    "\n",
    "axs[1].plot(sigmas,cents[0,:],'ks',markersize=8,markerfacecolor=(.3,.3,.3),label=r'$\\Delta$ to mean')\n",
    "axs[1].plot(sigmas,cents[1,:],'ko',markersize=8,markerfacecolor=(.8,.8,.8),label=r'$\\Delta$ to median')\n",
    "axs[1].legend(fontsize=12)\n",
    "axs[1].set(xlabel=r'$\\sigma$ parameter for exp(X$\\sigma$)',ylabel='Mean or median dist.')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Distance to H$_0$')\n",
    "\n",
    "axs[2].plot(sigmas,tvals[0,:],'ks',markersize=8,markerfacecolor=(.3,.3,.3),label='Parametric t')\n",
    "axs[2].plot(sigmas,tvals[1,:],'ko',markersize=8,markerfacecolor=(.8,.8,.8),label='Wilcoxon z')\n",
    "axs[2].legend(fontsize=12)\n",
    "axs[2].set(xlabel=r'$\\sigma$ parameter for exp(X$\\sigma$)',ylabel='z or t value')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Test stat. values')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex10a.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "njZ2jLWuvL-R"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# plot showing relationship between central tendency distances and test statistic values\n",
    "\n",
    "_,axs = plt.subplots(1,2,figsize=(8,3.5))\n",
    "\n",
    "axs[0].plot(cents[0,:],tvals[0,:],'ko',markersize=12,markerfacecolor=(.6,.6,.6))\n",
    "axs[0].set(xlabel='Distance: mean to .5',ylabel='T-value',xlim=[-.6,-.4],\n",
    "           title=r'$\\bf{A}$)  One-sample t-test')\n",
    "\n",
    "\n",
    "axs[1].plot(cents[1,:],tvals[1,:],'ko',markersize=12,markerfacecolor=(.6,.6,.6))\n",
    "axs[1].set(xlabel='Distance: median to .5',ylabel='Z-value',\n",
    "           title=r'$\\bf{B}$)  Wilcoxon test')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex10b.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Sqnr-HXGvL7g"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "mke93CTpvHMy"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 11"
   ],
   "metadata": {
    "id": "P5UlA-KKKgpj"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# params\n",
    "meanoffsets = np.linspace(0,2,71)\n",
    "samplesizes = np.arange(10,811,step=50)\n",
    "\n",
    "\n",
    "# initialize\n",
    "pvals = np.zeros((len(samplesizes),len(meanoffsets)))\n",
    "cohend = np.zeros((len(samplesizes),len(meanoffsets)))\n",
    "r2 = np.zeros((len(samplesizes),len(meanoffsets)))\n",
    "\n",
    "\n",
    "# loop over sample sizes\n",
    "for sidx,ssize in enumerate(samplesizes):\n",
    "\n",
    "  # loop over effect sizes\n",
    "  for eidx,effect in enumerate(meanoffsets):\n",
    "\n",
    "    # generate the data\n",
    "    X = np.random.normal(loc=effect,scale=1.5,size=(ssize))\n",
    "\n",
    "    # run the t-test and store the results\n",
    "    T = stats.ttest_1samp(X,0)\n",
    "    pvals[sidx,eidx] = T.pvalue\n",
    "\n",
    "    # Cohen's d\n",
    "    cohend[sidx,eidx] = np.abs(  np.mean(X)/np.std(X,ddof=1)  )\n",
    "\n",
    "    # R^2\n",
    "    r2[sidx,eidx] = T.statistic**2 / (T.statistic**2 + T.df)"
   ],
   "metadata": {
    "id": "QWPyL69JKglM"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "axs[0].plot(np.log(pvals+np.finfo(float).eps),cohend,'ko',markersize=8,markerfacecolor=(.8,.8,.8))\n",
    "axs[0].set(xlabel='log(p)',ylabel=\"Cohen's d\")\n",
    "axs[0].set_title(r\"$\\bf{A}$)  Cohen's d by p-values\")\n",
    "\n",
    "axs[1].plot(cohend,r2,'ko',markersize=8,markerfacecolor=(.8,.8,.8))\n",
    "axs[1].set(xlabel=\"Cohen's d\",ylabel=r'$R^2$')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Different effect size measures')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex11.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "-UNUjhVEKgf5"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "FBMBMPYdKgdU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 12"
   ],
   "metadata": {
    "id": "uXG7elbVvHOu"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# data reference\n",
    "\n",
    "# P. Cortez, A. Cerdeira, F. Almeida, T. Matos and J. Reis.\n",
    "# Modeling wine preferences by data mining from physicochemical properties. In Decision Support Systems, Elsevier, 47(4):547-553, 2009.\n",
    "\n",
    "# https://archive.ics.uci.edu/ml/datasets/Wine+Quality"
   ],
   "metadata": {
    "id": "zojHuxifujhH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "0ohXIxzt4_U2"
   },
   "source": [
    "url = \"https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv\"\n",
    "\n",
    "data = pd.read_csv(url,sep=';')\n",
    "data"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "WwRtOmUp5C80"
   },
   "source": [
    "# describe the data\n",
    "data.describe()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "RnOfpGIC5GOs"
   },
   "source": [
    "# list number of unique values per column\n",
    "for i in data.keys():\n",
    "  print(f'{i:>20} has {len(np.unique(data[i])):>3} unique values')"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "Ht-Xh9Db-ajg"
   },
   "source": [
    "# plot some data\n",
    "fig,ax = plt.subplots(1,figsize=(17,4))\n",
    "ax = sns.boxplot(data=data)\n",
    "ax.set_xticklabels(ax.get_xticklabels(),rotation=45)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# optionally remove outliers based on visual inspection (not used for the exercise)\n",
    "# data = data[data['total sulfur dioxide']<200]"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "BD3oTedE5Gfk"
   },
   "source": [
    "### z-score all variables except for quality\n",
    "\n",
    "# find the columns we want to normalize (all except quality)\n",
    "cols2zscore = data.keys()\n",
    "cols2zscore = cols2zscore.drop('quality')\n",
    "\n",
    "# make a copy of the dataframe to change\n",
    "dataz = data.copy()\n",
    "\n",
    "# z-score (written out for clarity)\n",
    "for col in cols2zscore:\n",
    "  meanval   = np.mean(data[col])\n",
    "  stdev     = np.std(data[col],ddof=1)\n",
    "  dataz[col] = (data[col]-meanval) / stdev\n",
    "\n",
    "# can also do more compactly\n",
    "#data[cols2zscore] = data[cols2zscore].apply(stats.zscore)\n",
    "\n",
    "dataz.describe()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "BRBQ7vdZ5GiU"
   },
   "source": [
    "# check the plot again\n",
    "fig,ax = plt.subplots(1,figsize=(15,4))\n",
    "ax = sns.boxplot(data=dataz)\n",
    "ax.set_xticklabels(ax.get_xticklabels(),rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# test each variable for normality\n",
    "\n",
    "# loop through all variables\n",
    "for col in cols2zscore:\n",
    "\n",
    "  # compute and print the test\n",
    "  Stest = stats.shapiro(dataz[col])\n",
    "  print(f'{col:>25}: p<{Stest.pvalue:.4f}')"
   ],
   "metadata": {
    "id": "1ZH9DAPoTjzs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualize the histograms\n",
    "_,axs = plt.subplots(3,4,figsize=(10,5))\n",
    "\n",
    "# loop through columns and create histograms\n",
    "for (idx,col),a in zip(enumerate(dataz.columns),axs.flatten()):\n",
    "  sns.histplot(data=dataz,x=col,ax=a,color=(.7,.7,.7))\n",
    "  a.set(xticks=[],yticks=[],ylabel='')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('ttest_ex12.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "cVYyk9gRTjtg"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "qjvI_6Su5Gk-"
   },
   "source": [
    "# binarize quality ratings into a new variable\n",
    "fig = plt.figure(figsize=(7,3))\n",
    "\n",
    "counts = dataz['quality'].value_counts()\n",
    "plt.bar(list(counts.keys()),counts,color=(.7,.7,.7),edgecolor='k')\n",
    "plt.xlabel('Quality rating')\n",
    "plt.ylabel('Count')\n",
    "plt.show()\n",
    "\n",
    "# create a new column for binarized (boolean) quality\n",
    "dataz['boolQuality'] = False\n",
    "dataz['boolQuality'][dataz['quality']>5] = True\n",
    "\n",
    "dataz[['quality','boolQuality']]"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "5Tqu6aRQAr9A"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 13"
   ],
   "metadata": {
    "id": "nzfttF3zT0Y4"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# run all t-tests and store the results\n",
    "\n",
    "# Note about this and the next code cell: You need to compute all p-values in order\n",
    "# to conduct the FDR correction. That's why I run the t-tests here and then report\n",
    "# the results in the following cell.\n",
    "\n",
    "\n",
    "\n",
    "# initialize results matrix as a dictionary\n",
    "results = {}\n",
    "\n",
    "\n",
    "# loop over column\n",
    "for col in cols2zscore:\n",
    "\n",
    "  # for convenience, extract the numerical variables\n",
    "  Xh = dataz[col][dataz['boolQuality']==True].values  # high rating\n",
    "  Xl = dataz[col][dataz['boolQuality']==False].values # low rating\n",
    "\n",
    "  # compute df\n",
    "  s1,s2 = np.var(Xh,ddof=1),np.var(Xl,ddof=1)\n",
    "  n1,n2 = len(Xh),len(Xl)\n",
    "  df_num = (s1/n1 + s2/n2)**2\n",
    "  df_den = s1**2/(n1**2*(n1-1)) + s2**2/(n2**2*(n2-1))\n",
    "\n",
    "  # run the t-test and store the results in a dictionary\n",
    "  tres = stats.ttest_ind(Xh,Xl,equal_var=False)\n",
    "  #tres = stats.mannwhitneyu(Xh,Xl) # uncomment for Mann-Whitney U test\n",
    "\n",
    "  results[col] = { 't' : tres.statistic,\n",
    "                   'p' : tres.pvalue,\n",
    "                   'df': df_num/df_den }\n"
   ],
   "metadata": {
    "id": "LxMb6PFJT0WN"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# need FDR correction function\n",
    "from statsmodels.stats.multitest import fdrcorrection\n",
    "\n",
    "# bonferroni threshold\n",
    "bonP = .05/len(cols2zscore)\n",
    "\n",
    "# FDR correction (don't need p-values, only keep outcome)\n",
    "fdrH = fdrcorrection([results[col]['p'] for col in cols2zscore],.05)[0]\n",
    "\n",
    "\n",
    "# loop through columns and report the results!\n",
    "for i,col in enumerate(cols2zscore):\n",
    "\n",
    "  # extract values\n",
    "  t  = results[col]['t']\n",
    "  p  = results[col]['p']\n",
    "  df = results[col]['df']\n",
    "\n",
    "  # determine if significant\n",
    "  issigB = [' ','*'][int(p<bonP)] # convert from bool to int to index\n",
    "  issigF = [' ','+'][int(fdrH[i])]\n",
    "\n",
    "  print(f'{col:>20}: t({df:.0f})={t:6.2f}, p={p:.4f}, {issigB}{issigF}')"
   ],
   "metadata": {
    "id": "3lshMiQJ2qHf"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "9eU97PQN2p_2"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}